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Abstract. This paper contains an analysis of a simple neural network that exhibits 
self-organized criticality. Such criticality follows from the combination of a simple 
neural network with an excitatory feedback loop that generates bistability, in 
combination with an anti-Hebbian synapse in its input pathway. Using the methods 
of statistical field theory, we show how one can formulate the stochastic dynamics 
of such a network as the action of a path integral, which we then investigate using 
renormalization group methods. The results indicate that the network exhibits 
hysteresis in switching back and forward between its two stable states, each of which 
loses its stability at a saddle-node bifurcation. The renormalization group analysis 
shows that the fluctuations in the neighborhood of such bifurcations have the signature 
of directed percolation. Thus the network states undergo the neural analog of a phase 
transition in the universality class of directed percolation. The network replicates 
precisely the behavior of the original sand-pile model of Bak, Tang & Wiesenfcld. 
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1. Introduction 



The idea of self- organized criticality (SOC) was introduced by Bak et al. (1988). Their 
paper immediately triggered an avalanche of papers on the topic, not the least of which 
was a connection with l/f- or scale-free noise. However it was not until another paper 
appeared, by Gil & Sornette (1996), which greatly clarified the dynamical prerequisites 
for achieving SOC, that a real understanding developed of the essential requirements 
for SOC: (1) an order-parameter equation for a dynamical system with a time-constant 
r Q , with stable states separated by a threshold, (2) a control-parameter equation with 
a time-constant r c , and (3) a steady driving force. In Bak ei.a/.'s classic example, the 
sand-pile model, the order parameter is the rate of flow of sand grains down a sand-pile, 
the control parameter is the sand-pile's slope, and the driving force is a steady flow of 
grains of sand onto the top of the pile. Gil and Sornette showed that if r G <C r c then the 
resulting avalanches of sand down the pile would have a scale-free distribution, whereas 
if t ^> r c then the distribution would also exhibit one or more large avalanches. 

In this paper we will analyze a neural network model which is in one-to-one 
correspondence with the Gil-Sornette SOC-model, and therefore also exhibits SOC. 



1.1. Neural network dynamics 

Consider first the mathematical representation of the dynamics of a neocortical slab 
comprising a single spatially homogeneous network of iV excitatory neurons. Such 
neurons make transitions from a quiescent state q to an activated state a at the rate a 
and back again to the quiescent state q at the rate a, as shown in figure 1. 




Figure 1. Neural state transitions 



Let P n (t) be the probability that a fraction n/N is active at time t. Then the 
probabilistic dynamics of such a slab can be formulated as a master equation of the 
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form: 

dP n (t) 



a[{n + l)P n+1 - nP n ] + (N - n + l)f[s(n - l)]P n _! 



- (N-n)/[s(n)]P n (1) 

where a is the rate at which activated neurons become quiescent, and s(n) is the total 
current or excitation driving each neuron in the population to fire at the rate a = f[s(n)]. 
We assume in this simplified model that each neuron receives a signal weighted by Wq/N 
from each of iV other neurons in the population, so that 

s(n) = Won + h (2) 

where h is an external current. 

Equation [I] can be extended to the spatially inhomogeneous case. Let n r /N r be 
the fraction of active cells at time t in the rth population of N r cells, and let P[n,i] be 
the probability of the configuration n = {n^, n 2 , • • • , n r , • • ■ , n^} existing at time t. The 
extended master equation then takes the form 

= a [{nr + l)P[n r+ , t] - n r P[n, t}} 

r 

+ ^[(N r - n r + l)f[s(n r - l)]P[n r .,t] 

r 

- {N r -n r )f[s(n r ))P[n,t)) (3) 

where n r ± = {n l5 n 2 , • • • , n r ± 1, • • • , uq} and there are a total of Q locally homogeneous 
populations. 

Following Van Kampen (1981) equation [3] can be rewritten using the one-step 
operators 

£?/(n) = f(n r ± 1) (4) 

i.e., 

= H £ r ~ l)nr + {e~ - l)(N r - n r )f[s(n r )}] P[n,t] (5) 

r 

Note that the total number of cells N r in the rth population comprises n r active 
cells and N r — n r = q r quiescent cells, so that 

n r + q r = N r (6) 

Thus equation |5] can be rewritten slightly in the form 

= E H S r ~ + (ST " VlrMrir)]] P[n,t] (7) 
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1.2. Annihilation and creation operators 

We now introduce Fock space annihilation and creation operators satisfying boson 
commutation rules 

[a r ,a\] = [q r ,qt] = 6(r - s) 
[a r ,a s ] = [a;,af s ] = 

[q r ,q a ] = [qt,qt]=0 (8) 
such that given neural vectors \n r ), \m r i) 

a\\n r ) = \n r + 1), a r \n r ) = n r \n r — 1) 

ql\m r ) = \m r + 1), g r |m r ) = m r \m r — 1) (9) 

where 

n r + m r = 1, \n r ) = (<f)?\0 r ), \m r ) = (gj)™|0 r ) (10) 
and |0 r ) is a fiducial or vacuum (empty) state such that 
al\0 r ) = \n r = l r ),a r |0 r ) = 

ql\0 r ) = \Tn r = l r ),g r |0 r ) = (11) 
and there exists a dual vector (0| such that 
(0 r |a r = (n r = l r \, (0 r \al = 0, 

(0 r \q r = (m r = l r |,(0 r |gj = 0, (12) 
It follows from equation [9] that 

a\.a r \n r ) = n r \n r ) 

i.e. n r is the eigenvalue of the operator a\.a ri which is therefore referred to as a number 
or number density operator. 

1.3. Bosonizing the master equation 

The master equation can now be bosonized by replacing the van Kampen one-step 
operators by bosonic equivalents, i.e. 

(S + - 1) = (qt - a\)a r 

{S- - 1) = {at - qt)q r (13) 
so eqn [7] becomes: 

2 = J2 W - 4)ar + (4 - ql)q r f[s(ala r )]] P[n, t] (14) 

r 

Let \P(t)) be a probability state vector satisfying 

|P(t)) = 5>M|n) (15) 
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Then equation [7] can be written in the form 

j t \p{t)) = H4 - 4H + (4 - <zJ:)<zr/K4a P )]] \p(t)) (ie) 

r 

or formally as 



where 



'P(f)> = -£|P(t)> (17) 



— H = ^ [a(4 — 4)°r + (4 — 4)?r/[s(4 a r)]] 
r 

is the quasi- Hamiltonian operator for the Markov process represented in equation [TJ 
i.^. From bosons to coherent states 



Equation [18] is a linear operator equation with formal solution 

\P(t))=exp[-H(t-t )\P{t )) 

We need to re-express this solution in terms of numbers rather than operators. This 
can be achieved by introducing coherent states. These were introduced by Schrodinger 
(1926) and first used extensively in coherent optics by Glauber (1963). We therefore 
introduce such states \<p r ) in the form 

\<p r ) = exp[--(p*(p r + <£v4]|0 r ) (19) 

where tp r is the right eigenvalue of a r , i.e. a r \cj) r ) = (p r \<f> r ). There is also a coherent 
state representation of q r in the form \9 r ) such that the right eigenvalue of q r is $ r , i.e. 
q r \6 r ) = d r \9 r ). In similar fashion (</v|4 = {(j) r \(p r where (p r , the complex conjugate of 
(p, is the left eigenvalue of aj, and similarly (Q T \ql = (0 r \$r, i>e. $ r is the left eigenvalue 
of q\. It follows that 

(cj) r \ala r \(j) r ) = (4> r \(f r ip r \(f) r } = ip r (f r (20) 

All this suggests that the operator quasi-Hamiltonian has a coherent state representation 
in the form 

~ Ti = [«($r - <fr)<Pr + (<^r ~ #r)$r/ [s((f r (p r )) (21) 

r 

Note in passing that operator products involving powers of the number operator 
a\a r must first be normal ordered, i.e. all creation operators a\ must preceed the 
annihilation operators a r , before coherent states can be introduced. For example the 
normal order form of exp[a^a] written as : exp[a^a] : is expanded as 

: expfa^a] : = 1 + (a)a) + —^{c^a) 2 + ■ ■ ■ 

= 1 + (a f a) + ^(a f a + a t2 a 2 ) H 

oo 1 I 
1=0 ' 3=0 
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where the Sij are Stirling numbers of the second kind. It follows that : expfa^a] : can 
be written as 

: expfa+a] := ^h k a ]k a k (23) 

A;=0 

where h k = Y,i «i,fe- 

The final preliminary step in this formulation is to take the continuum limit of the 



expression for H in equation 21, so that 

H I I d d xdt \a{d - <p)ip + {<p- ti)tff[s{0ip + <p)] (24) 




where (p r — > (p(pc,t) = (p etc., and the conjugate coherent state <p has been shifted to 

1 + 0. 

1.5. Dimensions and the density representation 

Before proceeding further we need to assign a dimension to each variable in equation [24} 
To do so we use a modified version of the convention used in particle physics so 
that [x] = L^ 1 , [t] = L~ 2 where L is the length scale used, whence [x 2 /t] = L°. 
{This generates a scaling found in Markov random walks and related processes such 
as stochastic neural activity.} Then [a] = L 2 , [if] = L d , [0\ = L°, [tpip] = L d , 
[f[s\] = [a] = L 2 . This last value of [f[s]] implies that the input current function 
s({pip + ip) = s(I) = kl where the constant k has the dimensions of inverse current. The 
net effect of such a choice leads to the required result that [H] = 0. 

To emphasize this choice we further transform the coherent-state quasi-Hamiltonian 
by introducing the density representation: 

(p — > exp[n] — 1, (p — > n[exp[—[n]] 

i? — > exp[p] — l,i? — > p[exp[— [p]] (25) 



so that equation 24 transforms into 



-H = J J d d xdt[a(exp(p-h) - l)n + (exp(n-p) - l)pf[s(n))) (26) 

Note that in the continuum limit the input current function s(n) = s(I) = kl = 
k(w~kn+h) where * is the spatial convolution operator, i.e. w-kn = J d d x'w('x— x')n(x', t). 

1.6. From the quasi-Hamiltonian to a neural Path Integral 

Using standard methods (Doi 1976, Peliti 1985) Buice & Cowan (2007) incorporated 
the quasi-Hamiltonian into the action of a Wiener path integral. This action takes the 
form: 

S(n) = J J d d xdt [nd t n + pd t p+ 

a(l — exp(— (n — p))n — (exp(n — p) — l)pf[s(n)}}(27) 
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The utility of this action is that it is part of the exponent of the moment generating 
functional for the statistical moments of the probability density P[n, t]. 
At an extremum 



5S(n) 



5S(n) 

n=0 S P 



= (28) 

p=0 



which leads to the mean-field equations 

d t n + an — pf[s(n)] = 0, d t p — an + pf[s(n)] = (29) 

whence 

n + p = p (30) 

where p is the (constant) packing density of excitatory neurons in the population. This 
is a mean- field result, so we replace p by p — n d in equation |27[ so as to generate 
the mean-field Wilson-Cowan equation (Wilson & Cowan 1973) for a single spatially 



organized population from the first variation of the action given in equation 27 and (in 



principle), all higher moments, and (finally) rewrite equation 27 in the form: 
S(n) = J J d d xdt [nd t n+ 

a(l - exp(-n))n - (exp(n) - l)(p - n d )f[s(n)}} (31) 

1.1. Renormalizing the path integral 



We expand n about its mean value (n) = n d , which satisfies equation 28 in the form: 
d t n d = -an d + (p - n d )f[s(n d )] (32) 

where 

s(n d ) = k(w -k n d + h d ) (33) 

Thus n — > n+n d , ft — >• n, since n d = 0. So s(n) —> s{n+n d ) and f[s{n)] —> f[s(n+n d )]. 
If follows that s{n + n d ) = k{w~k{n + n d ) + (h + h d )) = k(w-kn d + h d ) + k(w-kn + h)) = 
s{n d ) + s(n), and therefore f[s(n)) = f[s(n d ) + s(n)). We next expand f[s(n)] in 



a Taylor expansion about the mean-field value n d , noting that from equation A.l 

s(n) = k(w * n + h)) = k(Ln + h). In what immediately follows we assume that 
the external stimulus h(x,t) = 0. It follows that: 

f[s(n)\ = f[kL(n d + n)\ = f[kLn d ] + f^[kLn d ]kLn 

+ l -f^[kLn d ]{kLnf + --- (34) 



However because of normal ordering, equation [34] leads to the expression: 

m f(i) 

fl s ( n )] = ^2g m (kLn) m , where g m = ^ ~jj~ s ^m (35) 

m 1=0 

Since the leading terms of g m are proportional to f^ m \ and given the assumed form for 
f[s(ri)] to be such that p 1 ^ > and < 0, then g m > for m odd, and g m < for m 
even. 
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We also expand the functions exp(±n). The resulting action S(n) takes the form: 
S(n) = d d xdt [n(d t + a — (p — n c ])g\kV)n 



~y ( a + (P - n c i)gikL)n + n((p - n d )g 2 (kL) 2 )n 2 
n 2 

+ — {{p-n cl )g 2 {kL) 2 )n 2 + -- 



(36) 



It follows from the appendix that w 2 is small compared to Wo, so that in most 
expressions the terms proportional to V 2m n m can be neglected. However this is not 
always the case for m = 1. Thus the first term can be written approximately as 
h(d t + a- (p-n d )g 1 k(w + ^w 2 V 2 ))n = h(d t + p-DV 2 )n where p = a-(p-n d )gikwo 
and D = \{p — n c i)g\kw 2 . So the expression for the action is now reduced to the form: 

S(n) = J J d d xdt [h(d t + p- DV 2 )n - h 2 G x n 

+hG 2 n 2 + \h 2 G 2 n 2 + ■■■ (37) 

where Gi = l/2(a+(p— n c i)gikwo), and G 2 = (p—n c i)\g 2 \k 2 WQ. We need to demonstrate 
that the last term in S(n), i.e., ^n 2 G 2 n 2 , and all other terms, are irrelevant in the sense 
of the renormalization group. 

The renormalization group [RG] analysis is carried out via dimensional 
analysis. It can be shown that all the terms in S(n) are zero-dimensional when 
integrated over d- dimensional space and over time, i.e., [d d xdt] = L~( d+2 ) and 
[any term in the integrand] = L d+2 . However, as it stands [n] = L d , but [n] = L°, 
so that [hn] = L 0+d = L d . This is not suitable for the scaling analysis implemented in 
the RG process. We therefore introduce a new scaling, 

G\ _ / G 2 

Q-n, s=W— n (38) 

such that ss = hn where [G 2 /G\] = L~ d . The effect of this scaling is that both s and s 
have dimension L d / 2 . Let 



VG!G 2 = u, G 2 = 2t (39) 
The net effect of this scaling transformation is that 

S(s) = J J d d xdt [s(d t + p- DV 2 )s + uS(s - s)s + ts 2 s 2 + • • •] (40) 

But [t/u] = L~ d l 2 and therefore scales to zero as L — > oo under subsequent RG 

transforms. So asymptotically the terms rs 2 s 2 H become irrelevant in the RG sense. 

So finally 

S(s) = J J d d xdt[s{d t + p- DV 2 )s + u~s{s-s)s] (41) 
is the renormalized action of the large-scale neural activity of a single neural population. 
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1.8. Directed Percolation 

This action is well-known: it is called Reggeon Field Theory, and is found in directed 
percolation [DP] in random graphs, in contact processes, in high-energy nuclear physics, 
in bacterial colonies, all of which exhibit the characteristic properties of what is called 
a universality class, i.e., it is a phase transition with a universal scaling of important 
statistical exponents. It also shows up in branching and annihilating random walks, 
catalytic reactions, and interacting particles. Thus we have mapped the mathematics of 
large-scale neural activity in a single homogeneous neural population into a percolation 
problem in random graphs, or equivalently into a branching and annihilating random 
walk. A first version of this work was presented in Buice & Cowan (2007). A more 
extensive paper with many applications to neuroscience was presented in Buice & Cowan 
(2009). 

Here we note that there is an upper critical dimension at which directed percolation 
crosses over to mean-field behavior. This upper critical dimension is d — 4. What is 
the dimension of the neocortex? To answer this question we look at the number of 
synapses per neuron in the neocortex. Using estimates provided by Stevens (1989), 
this number is about 4 x 10 3 . Assuming the number of synapses in a terminal axonal 
arbor to be about 50, the number of neural neighbors per neuron is about 80, so the 
effective dimensionality of a neocortical hyper lattice is about d = 40. Thus the critical 
exponents characterizing the neural phase transition are the d = 4 exponents of directed 
percolation. These have been calculated by Abarbanel & Broznan (1974), Abarbanel 
et al. (1976), and Amati et al. (1976), and appear in the linear response of the neocortical 
model to an impulsive stimulus, known to mathematicians as the Green's function and 
to physicists as the propagator. This takes the form 



G(x-x',t-t') oc < 



(t-f)- 2 exp[-^^-fjL(t-f)\, fi>0 

(t-0- 2 exp[-^)l, n = (42) 



fi 2 Q \y/\iT\(t - 1>) 



4(t-f) _ 
\x — X 



fi<0 



where the cases fi > 0, /i = 0, and fi < correspond, respectively, to the sub-critical, 
critical, and super-critical propagators. They correspond, respectively, to solutions of 
the cable equation, the diffusion equation, and a nonlinear wave-equation. The critical 
propagator is thus the diffusion limit of Brownian motion. It turns out that there 
is a great deal of data supporting the hypothesis that the mean-field propagator of 
DP correctly describes the essential features of large-scale neocortical activity on many 
spatio-temporal scales. [See Burns (1951), Lampl et al. (1999), Nauhaus et al. (2009).] 
In addition data on the statistical structure of large-scale activity recorded in cortical 
slices by Beggs & Plenz (2003) supports the hypothesis. In particular, the avalanche-size 
distribution of spontaneous activity in cortical slices fits the DP hypothesis. 

The analysis can be extended to deal with a neural network comprising both 
excitatory and inhibitory neurons. However in this paper we describe how to incorporate 
synaptic plasticity into a network comprising only excitatory neurons, as a mechanism 
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that tunes the network so that it automatically reaches the critical point of the DP 
phase transition, thus exhibiting self- organized criticality. 



2. Incorporating synaptic plasticity 

Consider first a single excitatory population model with a fixed recurrent excitatory 
synapse we and an input H through an excitatory modifiable synapse wh- The mean- 




Figure 2. A recurrent excitatory network 



field neural equation for this is a version of equation 32, i.e. 
dnE 



dt 

where 



-a E n E + (1 - n E )o- E [se(Ie)] (43) 



se{Ie) = Ie/ Irh,e, I E = w E *n E + w H ~k n H (44) 
and the synapse % is modifiable with 

-gE{{n E - Pe,o - Pe,sWh) n H )t (45) 



duin 



dt 

where pe,o is a constant neural activity, pe,s is a constant (Vogels et al. 2011), and qe 
is a state-dependent rate function. 

This is also a mean-field equation in which the synaptic weight % is depressed by 
an anti-Hebbian mechanism, and potentiated by the input activity nu- The problem 
is to write an action that incorporates these equations. Note that the time scale of 
the growth and decay of neural activity is set by the constant ole-, whereas that of 
the growth and decay of synaptic plasticity is set by qe- Thus the ratio aE/gE is an 
important parameter. 

3. Deriving an action for synaptic plasticity 

The first problem is to develop an action for the modifiable synapse %. In order to do 



so we first note from equation A.2 that % scales with bn, the synaptic conductance or 



weight, and we can write an approximation to eqn. [45] in the form: 

= ~9e ( u e - Pe,o ~ pE,sk H b H ) n H (46) 
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We next reformulate the changes in bn as a Markov process with discrete states 
in continuous time. We therefore assume that bn is quantized in units A of synaptic 
weight, and similarly for &g. Thus 

b E = m E A, b H = m H A (47) 

and we look at the Markov process represented in figure |3j 



t + (m) t~(m+D 




Figure 3. Synaptic state transitions 

and we write a master equation for this process in the form 
dP(rriH, t) 



dt 



= r(m H + l)P(m H + l,t)- r(m H )P(m H , t) 
+ t + (m H - l)P(m H - 1, t) - t + (m H )P{m Hl t) 



(4£ 



where P(m,t) is the probability that the synaptic weight tuh =m at time t. and 

t + = 9e(Pe,o + pE,sk H fn H )n H 

t~ = g E n E n H (49) 
are the transition rates for the excitatory synapse m#. 

3.1. The van Kampen ladder operators 

We again introduce the van Kampen ladder operators 

E^m — m ± 1 
so that the master equation can be rewritten as: 

dP(mH, t) 



(50) 



dt 



[r(E+ H - 1) +t + (E~ H - 1)] P(m H ,t) 



(51) 



An examination of eqns. 48 50 indicates that only the transition rate t + contains a 
term proportional to m#. We will utilize this property in what follows. 
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3.2. From the Master Equation to the Action 

We now introduce bosonic annihilation and creation operators for m.#. Let such 
operators be denoted by s' and s respectively, and let \m) be a column vector 
representing the synaptic weight m such that 

s' \m) = \m + 1), s\m) = m\m — 1) (52) 

Note that s* and s satisfy the same commutation rules and equations that we introduced 
earlier, so that 

E+-l->8*- s t s = s t (l-s) 

E~ - 1 -> s - s j s = (1 - s f )s (53) 
So we consider the equation: 

d ^ = {E--l)\m) = {l-J)s\m) (54) 

We further note that (1 — s^)s is normal ordered. We therefore shift and s, 

s 1 " -> 1 + s, s -> s (55) 

where s and s are now coherent states, and introduce the density representation 

s — > exp(m) — 1, s— )-mexp(— m) (56) 

we find 

<9i|m) = — ss\m) 

= — (exp(m) — l)mexp(— m)\m) 

= — (1 — exp(— m))m\m) (57) 



Eqn. 21 tells us that the action S m must contain a term of the form 

-9EpE,sk H n H (^ - exp(-m H ))m H , 

a source term 

-9EPE,on H m H , 

and an interaction term of the form 

leading to an action of the form: 

S(m H ) = J J d d xdt [m H d t m H - g E pE,sk H n H {^ ~ exp(-m H ))m H 

-m H g E p Efi n H + rh H g E n E n H ] (5£ 



Using variational techniques, we can derive the mean-field equation [eqn 45 from the 
condition: 

SS(rh H ) 



SfriH 

For then we obtain the equation: 
drriH 



= (59) 



i.e., eqn. H6 



(lj ~9e (n E - p Efi ~ PE,skHm H ) n H (60) 
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3.3. Renormalizing the synaptic plasticity action 

We now proceed to renormalize the action S(ttih) just as we renormalized S(n E ). We 



therefore expand the exponential term in equation 58 and rewrite S{rriH) in the form: 

1 




S{m H ) = I I d d xdt 



rh H d t m H - H x m H m n n H + ^H^^itlhUh 



~H 2 rh H n H + m H g E n E n H } (61) 

where H x = g E p E)S ^u, and H 2 = QePe,o- 
We now introduce the scaling 

Sh = \ TT^H, s H = \ ~!r m H (62) 
V -"i V tl 2 

such that ShSh — rnHmH, and [H1/H2] = L~ d . This scaling is analogous to the scaling 

of n and n which we carried out earlier for neural activities. As before the effect of this 

scaling is that both sh and sh have dimension L d l 2 . 

Let 

VHiH 2 = u H , H x = 2r H (63) 



and recall that equation 38 scales n E to \/Gi/G 2 s E . 

Following the procedure outlined earlier we can calculate which terms in the 
transformed action S(sjj) become irrelevant under scaling transformations. The 
resulting renormalized synaptic plasticity action takes the form: 

S(s H ) = / / d d xdt[s H d t s H - UHS H n H } (64) 



4. Combining the actions 

It follows from this formulation that the full action for the coupled system of equations 
for the evolution of n E and m H can be obtained simply by adding the actions S{n E ) 
and S{t71h) together. The combined action therefore takes the form: 

S(n E ,rriH) = J J d d xdt [n E d t n E + a(l — exp(— n E ))n E 

-(exp(n E ) - l)(p - n Etd )f[s(n E )} + rh H d t m H 

-9EpE,sk H n H {^ - exp(-m H ))m H 

-m H g E p Efi n H + m H g E n E n H } (65) 

4-1. A simulation of the behavior of the combined mean-field equations 



The first variation of equation 65 generates the mean- field equations for n E and m#in 
the form: 

dn E 



dt 

drriff 
dt 



ot E n E + (1 - n E )a E [s E (I E )] 
9e {n E ~ Pe,o ~ pE,sk H m H ) n H (66) 
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where 

se(Ie) = k E (m E -kn E + m H -kn H ) (67) 
These equations can be simulated. The results are shown in figure 4. It will be seen that 
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Figure 4. Neural state transitions between a ground state and an excited state. 
Parameter values: m E = 3, nn = 3; a = 0.2. N* is the fixed-point value of he, 
and Wh is the magnitude of the anti-Hebbian synapse in the input path. 



in the "ground-state" of low values of N* = n* E the synaptic weight % (proportional to 
run), increases until it reaches the critical point at a saddle-node bifurcation, at which 
point N* becomes unstable and the system switches to the "excited-state". But then 
the anti-Hebbian term in the synaptic plasticity dynamics kicks in, and Wh declines 
until the excited-state fixed-point becomes unstable at the upper critical point, also 
at a saddle-node, and switches back to the ground-state fixed point, following which 
the hysteresis cycle starts over. This is a exact representation of the sand-pile model's 
behavior. This is another representation of SOC in a neural network. The reader should 
compare this with the synaptic mechanisms for achieving SOC described in Levina et al. 
(2007) and in Millman et al. (2010). 



4-2. Renormalizing the combined action 

To renormalize the combined action we follow the same procedure as before and 
expand the exponential functions exp[±n^], exp[±m#], and f[s{n E )}. Note that after 
normal ordering and collecting terms, /[«(%)] can be expanded in the extended form 
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f[s(n E )] = Ylm 9m(k E (L E n E + LhUh))" 1 - After scaling and dimensional analysis, the 
renormalized combined action takes the form: 

S(se,s h ) = J J d d xdt [s E (d t + he- D E V 2 )s E + u E ~s E {s E - ~s E )s E 

-v E s E n H + s H d t s H - u H s H n H ] (68) 

where u e ,Uh and v E are renormalized constants. It would appear that apart from 
the source term v E s E nn the coupled action is just the sum of the two uncoupled 
renormalized actions. This is indeed the case! All the addition terms which appear 
in the current function become irrelevant under renormalization, and do not effect the 
renormalized action representing n E . Similarly for m#. It follows that fluctuations 
in the activity n E in the neighborhood of the critical point /i E — 0, i.e. those in the 
fluctuation-driven regime, should be essentially those of directed percolation. 



5. Inhibitory synapses 

In case % is an inhibitory synapse the transitions t + (m) and t~(m) are reversed, so 
that: 

t + = g E n E n H 

t~ = 9e(pe,o + PE,skHm H )n H (69) 

Thus now only the transition rate t~ contains a term proportional to run- The effect of 
this is that the action for S m leads to the mean-field equation: 

= 9e (n E - p E ,o ~ PE,sk H m H ) n H (70) 

at an inhibitory synapse ran- 

Thus inhibitory feedforward synapses are Hebbian with stimulus dependent 
depression, whereas excitatory feedforward synapses are anti-Hebbian with stimulus 
dependent potentiation, and we have now formulated actions for feedforward excitatory 
and inhibitory synaptic plasticity, based on simple microscopic potentiation and 
depression processes, involving a unary variable, the synaptic weight m H . 



6. Conclusion 



In this paper we have indicated how one can formulate and analyze actions for a 
simple network of excitatory cells with an input coupled to the network via an activity- 
dependent modifiable synapse. This action allows, in principal, the computation of 
statistical moments of the fluctuating dynamics of the network. Previous work by Wilson 
& Cowan (1972) indicates that the mean-field dynamics of the network is bistable, and 
can generate a hysteresis loop. On coupling this dynamics to that of the modifiable 
feedforward synapse introduced, all the conditions for the achievement of SOC are 
present in the combined system. The simulation of the mean-field dynamics shows 
that SOC is indeed achieved. It only remains to simulate the behavior of the combined 
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system with intrinsic noise. Our prediction is that the fluctuations in the activity near 
the critical points of the system will exhibit the properties of directed percolation. This 
will be the subject of another paper. 
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Appendix 

Appendix A.l. Expanding the weighting function 
We approximate the convolution w * n as: 

= /^( X -x>(x<, ( ) 

- K + y W 2 V 2 + • • >(x', t) 

= Ln (A.l) 



w -kn 



where 



w(x) w(r) = b/a d e- r/a , a = r . (A.2) 

It follows that 

- - / = b rmh^ 

In case d — 3, wo — 8nb,W2 = 12w o- 2 , whence W2/1V0 = 12a 2 . This expansion of the 
weighting function w(x) is known as the moment expansion. 
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